1 📖 Introduction

This vignette demonstrates how to use CellChat (Jin, Plikus, and Nie 2025) to identify significant signaling changes across different biological conditions. This analysis is performed on single-cell RNA sequencing (scRNA-seq) data, and follows the general principles of scRNA-seq analysis workflows described in (Luecken and Theis 2019). The analysis can be performed using popular toolkits such as Seurat (Hao et al. 2024) and Scanpy (Wolf, Angerer, and Theis 2018). We will showcase CellChat’s diverse functionalities by applying it to scRNA-seq datasets from two distinct biological conditions: non-lesional (NL, normal) and lesional (LS, diseased) human skin from patients with atopic dermatitis. For a comprehensive comparison of different analysis methods, please refer to (Saelens et al. 2019)comparison.

Note: The two datasets in this tutorial have the same cell population compositions after joint clustering. For a guide on comparing datasets with different cell population compositions, please refer to the tutorial on Comparison analysis of multiple datasets with different cell type compositions.

CellChat utilizes a top-down approach to identify signaling changes at multiple levels. It begins with a high-level overview and progressively refines the analysis to reveal detailed insights into altered interactions, cell populations, signaling pathways, and ligand-receptor pairs.

2 📚 Load the required libraries

ptm = Sys.time()
library(CellChat)
library(patchwork)
library(ComplexHeatmap)
library(wordcloud)
library(ggplot2)
library(tidyverse)
library(ggalluvial)
library(ggrepel)
library(igraph)

2.1 📂 Create a directory to save figures

data.dir <- './comparison'
dir.create(data.dir)
# setwd(data.dir)

3 📦 Load CellChat object of each dataset and merge them together

First, run CellChat on each dataset individually. Then, merge the resulting CellChat objects. If you are using CellChat objects created with a version prior to 1.6.0, please update them using the updateCellChat function.

setwd("D:/Data_Analysis_Practice_2025/scRNAseq/CellChat-main")
cellchat.NL <- readRDS("data/cellchat_humanSkin_NL_analysed.rds")
cellchat.LS <- readRDS("data/cellchat_humanSkin_LS_analysed.rds")
object.list <- list(NL = cellchat.NL, LS = cellchat.LS)
cellchat <- mergeCellChat(object.list, add.names = names(object.list))
cellchat
execution.time = Sys.time() - ptm
print(as.numeric(execution.time, units = "secs"))
# Users can now export the merged CellChat object and the list of the two separate objects for later use
save(object.list, file = "cellchat_object.list_humanSkin_NL_LS.RData")
save(cellchat, file = "cellchat_merged_humanSkin_NL_LS.RData")

4 🔍 Identify altered interactions and cell populations

CellChat’s top-down analytical approach begins with a broad overview of the communication landscape and progressively delves into the finer details. This multi-level analysis identifies signaling changes in interactions, cell populations, pathways, and ligand-receptor pairs.

Initially, CellChat addresses the following key questions:

  • Is there an overall enhancement or suppression of cell-cell communication?
  • Which cell-cell interactions are significantly altered?
  • How do the major signaling sources and targets differ between conditions?

4.1 📊 Compare the total number of interactions and interaction strength

To determine if cell-cell communication is enhanced or suppressed, CellChat compares the total number of interactions and their strengths between the inferred communication networks of different biological conditions.

ptm = Sys.time()
gg1 <- compareInteractions(cellchat, show.legend = F, group = c(1,2))
gg2 <- compareInteractions(cellchat, show.legend = F, group = c(1,2), measure = "weight")

gg1 + gg2
compare interactions

Figure 1: compare Interactions

4.2 🧐 Compare the number of interactions and interaction strength among different cell populations

To pinpoint which cell populations exhibit significant changes in interaction, CellChat provides several visualization options for comparing interaction counts and strengths:

  • (A) Differential Circle Plot: Visualizes the differences in interactions between the two datasets.
  • (B) Differential Heatmap: Offers a detailed view of the interaction changes.
  • (C) Side-by-side Circle Plots: Compares the interaction networks for each dataset individually.
  • (D) Coarse-grained Analysis: Allows for the examination of differential interactions at a higher level by aggregating cell populations into defined groups.

4.2.1 ⭕ (A) Circle plot showing differential number of interactions or interaction strength among different cell populations across two datasets

The circle plot visualizes the differential number of interactions or their strengths between the two datasets. Edges colored in \(\color{red}{\text{red}}\) represent signaling that is increased in the second dataset, while \(\color{blue}{\text{blue}}\) edges indicate decreased signaling compared to the first dataset.

par(mfrow = c(1,2), xpd=TRUE)
netVisual_diffInteraction(cellchat, weight.scale = T)
netVisual_diffInteraction(cellchat, weight.scale = T, measure = "weight")

netVisual_diffInteraction your_second_image_alt_text

Figure 2: netVisual_diffInteraction

4.2.2 🔥 (B) Heatmap showing differential number of interactions or interaction strength among different cell populations across two datasets

For a more detailed view, CellChat can generate a heatmap of the differential number of interactions or their strengths. * The top bar plot shows the total incoming signaling changes for each cell population. * The right bar plot displays the total outgoing signaling changes. The height of the bars corresponds to the magnitude of the change between the two conditions. In the heatmap, \(\color{red}{\text{red}}\) indicates increased signaling in the second dataset, while \(\color{blue}{\text{blue}}\) represents decreased signaling.

gg1 <- netVisual_heatmap(cellchat)
gg2 <- netVisual_heatmap(cellchat, measure = "weight")

gg1 + gg2
netVisual_heatmap

Figure 3: netVisual_heatmap

4.2.3 🔵 (C) Circle plot showing the number of interactions or interaction strength among different cell populations across multiple datasets

The differential network analysis shown above is designed for pairwise comparisons. For studies involving more than two datasets, CellChat can visualize the number of interactions or their strengths for each dataset individually.

To ensure consistency in the visualizations across different datasets, CellChat normalizes the node sizes and edge weights by computing the maximum number of cells per cell group and the maximum interaction count (or weight) across all datasets.

weight.max <- getMaxWeight(object.list, attribute = c("idents","count"))
par(mfrow = c(1,2), xpd=TRUE)

for (i in 1:length(object.list)) {
  netVisual_circle(object.list[[i]]@net$count, weight.scale = T, label.edge= F, edge.weight.max = weight.max[2], edge.width.max = 12, title.name = paste0("Number of interactions - ", names(object.list)[i]))
}

netVisual_circle your_second_image_alt_text

Figure 4: netVisual_circle

4.2.4 🟣 (D) Circle plot showing the differential number of interactions or interaction strength among coarse cell types

To simplify the complex network and gain insights at a higher level, CellChat can aggregate cell-cell communication based on predefined cell groups. In this example, we categorize the cell populations into three main types and then re-merge the CellChat objects.

group.cellType <- c(rep("FIB", 4), rep("DC", 4), rep("TC", 4))
group.cellType <- factor(group.cellType, levels = c("FIB", "DC", "TC"))
object.list <- lapply(object.list, function(x) {mergeInteractions(x, group.cellType)})
cellchat <- mergeCellChat(object.list, add.names = names(object.list))

We can then visualize the number of interactions or their strengths between these aggregated cell types for each dataset.

weight.max <- getMaxWeight(object.list, slot.name = c("idents", "net", "net"), attribute = c("idents","count", "count.merged"))

par(mfrow = c(1,2), xpd=TRUE)
for (i in 1:length(object.list)) {
  netVisual_circle(object.list[[i]]@net$count.merged, weight.scale = T, label.edge= T, edge.weight.max = weight.max[3], edge.width.max = 12, title.name = paste0("Number of interactions - ", names(object.list)[i]))
}

netVisual_circle_merged your_second_image_alt_text

Figure 5: netVisual_circle_merged

Similarly, a circle plot can be used to show the differential number of interactions or their strengths between the aggregated cell types. Red edges indicate increased signaling in the second dataset, while blue edges represent decreased signaling.

par(mfrow = c(1,2), xpd=TRUE)
netVisual_diffInteraction(cellchat, weight.scale = T, measure = "count.merged", label.edge = T)
netVisual_diffInteraction(cellchat, weight.scale = T, measure = "weight.merged", label.edge = T)

netVisual_circle_merged your_second_image_alt_text

Figure 6: netVisual_diffInteraction_merged

4.3 🎯 Compare the major sources and targets in a 2D space

By plotting the outgoing versus incoming interaction strengths in a 2D space, we can readily identify cell populations that have undergone significant changes in their signaling roles (as senders or receivers) between the datasets.

Follow option (A) to identify the cell populations with the most significant changes, and option (B) to inspect the signaling changes for specific cell populations of interest.

4.3.1 (A) Identify cell populations with significant changes in sending or receiving signals

num.link <- sapply(object.list, function(x) {rowSums(x@net$count) + colSums(x@net$count)-diag(x@net$count)})
weight.MinMax <- c(min(num.link), max(num.link)) # control the dot size in the different datasets
gg <- list()
for (i in 1:length(object.list)) {
  object.list[[i]] <- netAnalysis_computeCentrality(object.list[[i]])
  gg[[i]] <- netAnalysis_signalingRole_scatter(object.list[[i]], title = names(object.list)[i], weight.MinMax = weight.MinMax)
}

patchwork::wrap_plots(plots = gg)
netAnalysis_signalingRole_scatter

Figure 7: netAnalysis_signalingRole_scatter

The scatter plot reveals that Inflam.DC and cDC1 emerge as major signaling sources and targets in the lesional (LS) condition compared to the non-lesional (NL) one. Additionally, fibroblast populations become more prominent as signaling sources in the LS condition.

4.3.2 ⚡️ (B) Identify the signaling changes of specific cell populations

We can now delve deeper and identify the specific signaling pathways that are altered for Inflam.DC and cDC1 between the NL and LS conditions.

# compute centrality for each dataset and merge them
object.list <- lapply(object.list, function(x) {
  x <- netAnalysis_computeCentrality(x)
  return(x)
})
cellchat <- mergeCellChat(object.list, add.names = names(object.list))
gg1 <- netAnalysis_signalingChanges_scatter(cellchat, idents.use = "Inflam. DC", signaling.exclude = "MIF")
gg2 <- netAnalysis_signalingChanges_scatter(cellchat, idents.use = "cDC1", signaling.exclude = c("MIF"))

patchwork::wrap_plots(plots = list(gg1,gg2))

execution.time = Sys.time() - ptm
print(as.numeric(execution.time, units = "secs"))
netAnalysis_signalingChanges_scatter

Figure 8: netAnalysis_signalingChanges_scatter

5 🧠 Identify altered signaling with distinct network architecture and interaction strength

CellChat can identify signaling pathways with significant alterations in their network architecture or interaction strength. It can also classify signaling groups based on their functional or structural similarity across multiple biological conditions.

5.1 🕸️ Identify signaling networks with larger (or less) difference as well as signaling groups based on their functional/structure similarity

CellChat performs joint manifold learning and classification of the inferred communication networks based on their functional and topological similarity across different conditions. This analysis is applicable to more than two datasets.

By quantifying the similarity between the cellular communication networks of signaling pathways across conditions, this analysis highlights potentially altered signaling pathways. CellChat adopts the concept of network rewiring from network biology and hypothesizes that the difference between communication networks may affect biological processes across conditions. UMAP is used for visualizing signaling relationships and interpreting the outputs in an intuitive way.

Functional similarity: A high degree of functional similarity indicates that major senders and receivers are similar, suggesting that the two signaling pathways or ligand-receptor pairs exhibit similar or redundant roles. NB: Functional similarity analysis is not applicable to multiple datasets with different cell type compositions.

Structural similarity: Structural similarity is used to compare signaling network structures, without considering the similarity of senders and receivers. NB: Structural similarity analysis is applicable to multiple datasets with either the same or vastly different cell type compositions.

Here, we can run the manifold and classification learning analysis based on functional similarity because the two datasets have the same cell type composition.

5.1.1 ⚙️ Identify signaling groups based on their functional similarity

ptm = Sys.time()

cellchat <- computeNetSimilarityPairwise(cellchat, type = "functional")
cellchat <- netEmbedding(cellchat, type = "functional")
cellchat <- netClustering(cellchat, type = "functional")
# Visualization in 2D-space

netVisual_embeddingPairwise(cellchat, type = "functional", label.size = 3.5)

5.1.2 🏗️ Identify signaling groups based on structure similarity

cellchat <- computeNetSimilarityPairwise(cellchat, type = "structural")
cellchat <- netEmbedding(cellchat, type = "structural")
cellchat <- netClustering(cellchat, type = "structural")
# Visualization in 2D-space

netVisual_embeddingPairwise(cellchat, type = "structural", label.size = 3.5)
netVisual_embeddingPairwiseZoomIn(cellchat, type = "structural", nCol = 2)
netVisual_embeddingPairwise

Figure 9: netVisual_embeddingPairwise

5.1.3 📏 Compute and visualize the pathway distance in the learned joint manifold

CellChat can identify signaling networks with larger or smaller differences based on their Euclidean distance in the shared two-dimensional space. A larger distance implies a greater difference between the communication networks of the two datasets in terms of either functional or structural similarity. Note that we only compute the distance for signaling pathways that are present in both datasets. Signaling pathways identified in only one dataset are not considered. For comparisons involving more than three datasets, you can perform pairwise comparisons by modifying the comparison parameter in the rankSimilarity function.

rankSimilarity(cellchat, type = "functional")
rankSimilarity

Figure 10: rankSimilarity

5.2 🔄 Identify altered signaling with distinct interaction strength

By comparing the information flow (interaction strength) of each signaling pathway, CellChat identifies pathways that are turned off, decreased, turned on, or increased in one condition compared to another. You can identify altered signaling pathways or ligand-receptor pairs based on the overall information flow (Option A) or by examining outgoing/incoming signaling patterns (Option B).

5.2.1 🌊 (A) Compare the overall information flow of each signaling pathway or ligand-receptor pair

CellChat can identify conserved and context-specific signaling pathways by comparing the information flow for each pathway. The information flow is defined as the sum of communication probabilities among all pairs of cell groups in the inferred network (i.e., the total weights in the network).

This bar chart can be plotted in a stacked or non-stacked mode. Significant signaling pathways are ranked based on the differences in their overall information flow between the NL and LS skin conditions. When do.stat = TRUE, a paired Wilcoxon test is performed to assess the significance of the difference in information flow between the two conditions. The top signaling pathways colored red are enriched in NL skin, and those colored green are enriched in LS skin.

gg1 <- rankNet(cellchat, mode = "comparison", measure = "weight", sources.use = NULL, targets.use = NULL, stacked = T, do.stat = TRUE)
gg2 <- rankNet(cellchat, mode = "comparison", measure = "weight", sources.use = NULL, targets.use = NULL, stacked = F, do.stat = TRUE)

gg1 + gg2
rankNet

Figure 11: rankNet

5.2.2 📤📥 (B) Compare outgoing (or incoming) signaling patterns associated with each cell population

The previous analysis summarized the information from both outgoing and incoming signaling. CellChat can also compare the outgoing or incoming signaling patterns between two datasets, allowing for the identification of signaling pathways or ligand-receptors that exhibit different signaling patterns. We can combine all identified signaling pathways from the different datasets and compare them side-by-side, including outgoing, incoming, and overall signaling (aggregated from outgoing and incoming).

CellChat uses a heatmap to show the contribution of signals (signaling pathways or ligand-receptor pairs) to cell groups in terms of outgoing or incoming signaling. In this heatmap, the colorbar represents the relative signaling strength of a pathway across cell groups (values are row-scaled). The top colored bar plot shows the total signaling strength of a cell group by summarizing all pathways in the heatmap. The right grey bar plot shows the total signaling strength of a pathway by summarizing all cell groups.

i = 1
# combining all the identified signaling pathways from different datasets 
pathway.union <- union(object.list[[i]]@netP$pathways, object.list[[i+1]]@netP$pathways)
netAnalysis_computeCentrality(object.list[[i]])
ht1 = netAnalysis_signalingRole_heatmap(object.list[[i]], pattern = "outgoing", signaling = pathway.union, title = names(object.list)[i], width = 8, height = 8)
ht2 = netAnalysis_signalingRole_heatmap(object.list[[i+1]], pattern = "outgoing", signaling = pathway.union, title = names(object.list)[i+1], width = 8, height = 8)

draw(ht1 + ht2, ht_gap = unit(0.5, "cm"))

execution.time = Sys.time() - ptm
print(as.numeric(execution.time, units = "secs"))
complexheatmap

Figure 12: complexheatmap

ht1 = netAnalysis_signalingRole_heatmap(object.list[[i]], pattern = "incoming", signaling = pathway.union, title = names(object.list)[i], width = 8, height = 8, color.heatmap = "GnBu")
ht2 = netAnalysis_signalingRole_heatmap(object.list[[i+1]], pattern = "incoming", signaling = pathway.union, title = names(object.list)[i+1], width = 8, height = 8, color.heatmap = "GnBu")

draw(ht1 + ht2, ht_gap = unit(0.5, "cm"))
complexheatmap_incoming

Figure 13: complexheatmap_incoming

ht1 = netAnalysis_signalingRole_heatmap(object.list[[i]], pattern = "all", signaling = pathway.union, title = names(object.list)[i], width = 8, height = 8, color.heatmap = "OrRd")
ht2 = netAnalysis_signalingRole_heatmap(object.list[[i+1]], pattern = "all", signaling = pathway.union, title = names(object.list)[i+1], width = 8, height = 8, color.heatmap = "OrRd")

draw(ht1 + ht2, ht_gap = unit(0.5, "cm"))
complexheatmap_all

Figure 14: complexheatmap_all

6 💡 Identify the up-gulated and down-regulated signaling ligand-receptor pairs

6.1 🎲 Identify dysfunctional signaling by comparing the communication probabities

CellChat can compare the communication probabilities mediated by L-R pairs from certain cell groups to other cell groups. This can be done by setting comparison in the function netVisual_bubble.

ptm = Sys.time()

netVisual_bubble(cellchat, sources.use = 4, targets.use = c(5:11),  comparison = c(1, 2), angle.x = 45)
bubble_plot

Figure 15: bubble_plot

Moreover, CellChat can identify the up-regulated (increased) and down-regulated (decreased) signaling ligand-receptor pairs in one dataset compared to the other dataset. This can be done by specifying max.dataset and min.dataset in the function netVisual_bubble. The increased signaling means these signaling have higher communication probability (strength) in the second dataset compared to the first dataset. The ligand-receptor pairs shown in the bubble plot can be accessed via gg1$data.

gg1 <- netVisual_bubble(cellchat, sources.use = 4, targets.use = c(5:11),  comparison = c(1, 2), max.dataset = 2, title.name = "Increased signaling in LS", angle.x = 45, remove.isolate = T)
gg2 <- netVisual_bubble(cellchat, sources.use = 4, targets.use = c(5:11),  comparison = c(1, 2), max.dataset = 1, title.name = "Decreased signaling in LS", angle.x = 45, remove.isolate = T)

gg1 + gg2

NB: The ligand-receptor pairs shown in the bubble plot can be accessed via signaling.LSIncreased = gg1$data.

bubble_plot_up_down

Figure 16: bubble_plot_up_down

6.2 🔬 Identify dysfunctional signaling by using differential expression analysis

The above method for identifying the upgulated and down-regulated signaling is perfomed by comparing the communication probability between two datasets for each L-R pair and each pair of cell groups. Alternative, we can identify the upgulated and down-regulated signaling ligand-receptor pairs based on the differential expression analysis (DEA). Specifically, we perform differential expression analysis between two biological conditions (i.e., NL and LS) for each cell group, and then obtain the upgulated and down-regulated signaling based on the fold change of ligands in the sender cells and receptors in the receiver cells.

Of note, users may observe the same LR pairs appearing in both the up-regulated and down-regulated results due to the fact that DEA between conditions is performed for each cell group. To perform DEA between conditions by ignoring cell group information, users can set group.DE.combined = TRUE in identifyOverExpressedGenes for CellChat v2.1.1.

# define a positive dataset, i.e., the dataset with positive fold change against the other dataset
pos.dataset = "LS"
# define a char name used for storing the results of differential expression analysis
features.name = paste0(pos.dataset, ".merged")

# perform differential expression analysis 
# Of note, compared to CellChat version < v2, CellChat v2 now performs an ultra-fast Wilcoxon test using the presto package, which gives smaller values of logFC. Thus we here set a smaller value of thresh.fc compared to the original one (thresh.fc = 0.1). Users can also provide a vector and dataframe of customized DEGs by modifying the cellchat@var.features$LS.merged and cellchat@var.features$LS.merged.info. 

cellchat <- identifyOverExpressedGenes(cellchat, group.dataset = "datasets", pos.dataset = pos.dataset, features.name = features.name, only.pos = FALSE, thresh.pc = 0.1, thresh.fc = 0.05,thresh.p = 0.05, group.DE.combined = FALSE) 

# map the results of differential expression analysis onto the inferred cell-cell communications to easily manage/subset the ligand-receptor pairs of interest
net <- netMappingDEG(cellchat, features.name = features.name, variable.all = TRUE)
# extract the ligand-receptor pairs with upregulated ligands in LS
net.up <- subsetCommunication(cellchat, net = net, datasets = "LS",ligand.logFC = 0.05, receptor.logFC = NULL)
# extract the ligand-receptor pairs with upregulated ligands and upregulated receptors in NL, i.e.,downregulated in LS
net.down <- subsetCommunication(cellchat, net = net, datasets = "NL",ligand.logFC = -0.05, receptor.logFC = NULL)

Since the signaling genes in the net.up and net.down might be complex with multi-subunits, we can do further deconvolution to obtain the individual signaling genes.

gene.up <- extractGeneSubsetFromPair(net.up, cellchat)
gene.down <- extractGeneSubsetFromPair(net.down, cellchat)

Users can also find all the significant outgoing/incoming/both signaling according to the customized features and cell groups of interest

df <- findEnrichedSignaling(object.list[[2]], features = c("CCL19", "CXCL12"), idents = c("Inflam. FIB", "COL11A1+ FIB"), pattern ="outgoing")

6.3 🎨 Visualize the identified up-regulated and down-regulated signaling ligand-receptor pairs

CellChat can visualize the identified up-regulated and down-regulated signaling ligand-receptor pairs using bubble plot (option A), chord diagram (option B) or wordcloud (option C).

6.3.1 🫧 (A) Bubble plot

We then visualize the upgulated and down-regulated signaling ligand-receptor pairs using bubble plot or chord diagram.

pairLR.use.up = net.up[, "interaction_name", drop = F]
gg1 <- netVisual_bubble(cellchat, pairLR.use = pairLR.use.up, sources.use = 4, targets.use = c(5:11), comparison = c(1, 2),  angle.x = 90, remove.isolate = T,title.name = paste0("Up-regulated signaling in ", names(object.list)[2]))
pairLR.use.down = net.down[, "interaction_name", drop = F]
gg2 <- netVisual_bubble(cellchat, pairLR.use = pairLR.use.down, sources.use = 4, targets.use = c(5:11), comparison = c(1, 2),  angle.x = 90, remove.isolate = T,title.name = paste0("Down-regulated signaling in ", names(object.list)[2]))

gg1 + gg2
bubble_plot_up_down_DEA

Figure 17: bubble_plot_up_down_DEA

6.3.2 🎶 (B) Chord diagram

Visualize the upgulated and down-regulated signaling ligand-receptor pairs using Chord diagram

# Chord diagram

par(mfrow = c(1,2), xpd=TRUE)
netVisual_chord_gene(object.list[[2]], sources.use = 4, targets.use = c(5:11), slot.name = 'net', net = net.up, lab.cex = 0.8, small.gap = 3.5, title.name = paste0("Up-regulated signaling in ", names(object.list)[2]))
netVisual_chord_gene(object.list[[1]], sources.use = 4, targets.use = c(5:11), slot.name = 'net', net = net.down, lab.cex = 0.8, small.gap = 3.5, title.name = paste0("Down-regulated signaling in ", names(object.list)[2]))

6.3.3 ☁️ (C) Wordcloud plot

Visualize the enriched ligands, signaling,or ligand-receptor pairs in one condition compared to another condition using wordcloud

# visualize the enriched ligands in the first condition
par(mfrow = c(1,2), xpd=TRUE)
computeEnrichmentScore(net.down, species = 'human', variable.both = TRUE)
# visualize the enriched ligands in the second condition
computeEnrichmentScore(net.up, species = 'human', variable.both = TRUE)
wordcloud_up_down_DEA

Figure 18: wordcloud_up_down_DEA

7 👀 Visually compare cell-cell communication using Hierarchy plot, Circle plot or Chord diagram

Similar to the CellChat analysis of individual dataset, CellChat can visually compare cell-cell communication networks using hierarchy plot, circle plot, chord diagram, or heatmap. More details on the visualization can be found in the CellChat analysis of individual dataset.

Edge color/weight, node color/size/shape: In all visualization plots, edge colors are consistent with the sources as sender, and edge weights are proportional to the interaction strength. Thicker edge line indicates a stronger signal. In the Hierarchy plot and Circle plot, circle sizes are proportional to the number of cells in each cell group. In the hierarchy plot, solid and open circles represent source and target, respectively. In the Chord diagram, the inner thinner bar colors represent the targets that receive signal from the corresponding outer bar. The inner bar size is proportional to the signal strength received by the targets. Such inner bar is helpful for interpreting the complex chord diagram. Note that there exist some inner bars without any chord for some cell groups, please just igore it because this is an issue that has not been addressed by circlize package.

pathways.show <- c("CXCL") 
weight.max <- getMaxWeight(object.list, slot.name = c("netP"), attribute = pathways.show) # control the edge weights across different datasets

par(mfrow = c(1,2), xpd=TRUE)
for (i in 1:length(object.list)) {
  netVisual_aggregate(object.list[[i]], signaling = pathways.show, layout = "circle", edge.weight.max = weight.max[1], edge.width.max = 10, signaling.name = paste(pathways.show, names(object.list)[i]))
}
hierarchy_CXCL

Figure 19: hierarchy_CXCL

pathways.show <- c("CXCL") 

par(mfrow = c(1,2), xpd=TRUE)
ht <- list()
for (i in 1:length(object.list)) {
  ht[[i]] <- netVisual_heatmap(object.list[[i]], signaling = pathways.show, color.heatmap = "Reds",title.name = paste(pathways.show, "signaling ",names(object.list)[i]))
}
ComplexHeatmap::draw(ht[[1]] + ht[[2]], ht_gap = unit(0.5, "cm"))
heatmap_CXCL

Figure 20: heatmap_CXCL

# Chord diagram
pathways.show <- c("CXCL") 

par(mfrow = c(1,2), xpd=TRUE)
for (i in 1:length(object.list)) {
  netVisual_aggregate(object.list[[i]], signaling = pathways.show, layout = "chord", signaling.name = paste(pathways.show, names(object.list)[i]))
}
chord_CXCL

Figure 21: chord_CXCL

For the chord diagram, CellChat has an independent function netVisual_chord_cell to flexibly visualize the signaling network by adjusting different parameters in the circlize package. For example, we can define a named char vector group to create multiple-group chord diagram, e.g., grouping cell clusters into different cell types.

# Chord diagram
group.cellType <- c(rep("FIB", 4), rep("DC", 4), rep("TC", 4)) # grouping cell clusters into fibroblast, DC and TC cells
names(group.cellType) <- levels(object.list[[1]]@idents)
pathways.show <- c("CXCL") 

par(mfrow = c(1,2), xpd=TRUE)
for (i in 1:length(object.list)) {
  netVisual_chord_cell(object.list[[i]], signaling = pathways.show, group = group.cellType, title.name = paste0(pathways.show, " signaling network - ", names(object.list)[i]))
}
 chord_cellType_CXCL

Figure 22: chord_cellType_CXCL

Using chord diagram, CellChat provides two functions netVisual_chord_cell and netVisual_chord_gene for visualizing cell-cell communication with different purposes and different levels. netVisual_chord_cell is used for visualizing the cell-cell communication between different cell groups (where each sector in the chord diagram is a cell group), and netVisual_chord_gene is used for visualizing the cell-cell communication mediated by mutiple ligand-receptors or signaling pathways (where each sector in the chord diagram is a ligand, receptor or signaling pathway.)

par(mfrow = c(1, 2), xpd=TRUE)
# compare all the interactions sending from Inflam.FIB to DC cells
for (i in 1:length(object.list)) {
  netVisual_chord_gene(object.list[[i]], sources.use = 4, targets.use = c(5:8), lab.cex = 0.5, title.name = paste0("Signaling from Inflam.FIB - ", names(object.list)[i]))
}
chord_gene

Figure 23: chord_gene

# compare all the interactions sending from fibroblast to inflamatory immune cells
par(mfrow = c(1, 2), xpd=TRUE)
for (i in 1:length(object.list)) {
  netVisual_chord_gene(object.list[[i]], sources.use = c(1,2, 3, 4), targets.use = c(8,10),  title.name = paste0("Signaling received by Inflam.DC and .TC - ", names(object.list)[i]), legend.pos.x = 10)
}
chord_gene_fibroblast_inflam

Figure 24: chord_gene_fibroblast_inflam

# show all the significant signaling pathways from fibroblast to immune cells
par(mfrow = c(1, 2), xpd=TRUE)
for (i in 1:length(object.list)) {
  netVisual_chord_gene(object.list[[i]], sources.use = c(1,2,3,4), targets.use = c(5:11),slot.name = "netP", title.name = paste0("Signaling pathways sending from fibroblast - ", names(object.list)[i]), legend.pos.x = 10)
}

NB: Please ignore the note when generating the plot such as “Note: The first link end is drawn out of sector ‘MIF’.”. If the gene names are overlapped, you can adjust the argument small.gap by decreasing the value.

chord_gene_fibroblast_immune

Figure 25: chord_gene_fibroblast_immune

8 🧬 Compare the signaling gene expression distribution between different datasets

We can plot the gene expression distribution of signaling genes related to L-R pairs or signaling pathway using a Seurat wrapper function plotGeneExpression.

cellchat@meta$datasets = factor(cellchat@meta$datasets, levels = c("NL", "LS")) # set factor level

plotGeneExpression(cellchat, signaling = "CXCL", split.by = "datasets", colors.ggplot = T, type = "violin")

execution.time = Sys.time() - ptm
print(as.numeric(execution.time, units = "secs"))
plotGeneExpression_CXCL

Figure 26: plotGeneExpression_CXCL

9 ℹ️ Session Info

sessionInfo()
#> R version 4.5.0 (2025-04-11 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows 11 x64 (build 26100)
#> 
#> Matrix products: default
#>   LAPACK version 3.12.1
#> 
#> locale:
#> [1] LC_COLLATE=English_United States.utf8  LC_CTYPE=English_United States.utf8    LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C                           LC_TIME=English_United States.utf8    
#> 
#> time zone: America/New_York
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] grid      stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] ggrepel_0.9.6         ggalluvial_0.12.5     lubridate_1.9.4       forcats_1.0.1         stringr_1.5.2        
#>  [6] purrr_1.1.0           readr_2.1.5           tidyr_1.3.1           tibble_3.3.0          tidyverse_2.0.0      
#> [11] wordcloud_2.6         RColorBrewer_1.1-3    ComplexHeatmap_2.24.1 patchwork_1.3.2       CellChat_2.2.0       
#> [16] Biobase_2.68.0        BiocGenerics_0.54.0   generics_0.1.4        ggplot2_4.0.0         igraph_2.1.4         
#> [21] dplyr_1.1.4           htmltools_0.5.8.1     kableExtra_1.4.0      knitr_1.50           
#> 
#> loaded via a namespace (and not attached):
#>   [1] pbapply_1.7-4         rlang_1.1.6           magrittr_2.0.4        clue_0.3-66           GetoptLong_1.0.5     
#>   [6] gridBase_0.4-7        matrixStats_1.5.0     compiler_4.5.0        png_0.1-8             systemfonts_1.3.1    
#>  [11] vctrs_0.6.5           reshape2_1.4.4        sysfonts_0.8.9        pkgconfig_2.0.3       shape_1.4.6.1        
#>  [16] crayon_1.5.3          fastmap_1.2.0         backports_1.5.0       promises_1.3.3        rmarkdown_2.30       
#>  [21] tzdb_0.5.0            network_1.19.0        xfun_0.53             showtext_0.9-7        cachem_1.1.0         
#>  [26] jsonlite_2.0.0        later_1.4.4           irlba_2.3.5.1         broom_1.0.10          parallel_4.5.0       
#>  [31] cluster_2.1.8.1       R6_2.6.1              bslib_0.9.0           stringi_1.8.7         reticulate_1.43.0    
#>  [36] car_3.1-3             parallelly_1.45.1     jquerylib_0.1.4       Rcpp_1.1.0            assertthat_0.2.1     
#>  [41] iterators_1.0.14      future.apply_1.20.0   klippy_0.0.0.9500     IRanges_2.42.0        FNN_1.1.4.1          
#>  [46] timechange_0.3.0      httpuv_1.6.16         Matrix_1.7-4          tidyselect_1.2.1      abind_1.4-8          
#>  [51] rstudioapi_0.17.1     dichromat_2.0-0.1     yaml_2.3.10           doParallel_1.0.17     codetools_0.2-20     
#>  [56] listenv_0.9.1         lattice_0.22-7        plyr_1.8.9            shiny_1.11.1          withr_3.0.2          
#>  [61] S7_0.2.0              coda_0.19-4.1         evaluate_1.0.5        future_1.67.0         xml2_1.4.0           
#>  [66] circlize_0.4.16       ggpubr_0.6.1          pillar_1.11.1         BiocManager_1.30.26   carData_3.0-5        
#>  [71] rngtools_1.5.2        foreach_1.5.2         stats4_4.5.0          hms_1.1.3             S4Vectors_0.46.0     
#>  [76] scales_1.4.0          NMF_0.28              ggnetwork_0.5.14      globals_0.18.0        xtable_1.8-4         
#>  [81] glue_1.8.0            tools_4.5.0           BiocNeighbors_2.2.0   RSpectra_0.16-2       ggsignif_0.6.4       
#>  [86] registry_0.5-1        cowplot_1.2.0         colorspace_2.1-2      showtextdb_3.0        Formula_1.2-5        
#>  [91] cli_3.6.5             textshaping_1.0.3     viridisLite_0.4.2     svglite_2.2.1         gtable_0.3.6         
#>  [96] rstatix_0.7.2         sass_0.4.10           digest_0.6.37         sna_2.8               rjson_0.2.23         
#> [101] farver_2.1.2          lifecycle_1.0.4       statnet.common_4.12.0 GlobalOptions_0.1.2   mime_0.13

🔗 References

Hao, Y., T. Stuart, M. H. Kowalski, S. Choudhary, P. Hoffman, A. Hartman, A. Srivastava, et al. 2024. “Dictionary Learning for Integrative, Multimodal and Scalable Single-Cell Analysis.” Journal Article. Nat Biotechnol 42 (2): 293–304. https://doi.org/10.1038/s41587-023-01767-y.
Jin, S., M. V. Plikus, and Q. Nie. 2025. “CellChat for Systematic Analysis of Cell-Cell Communication from Single-Cell Transcriptomics.” Journal Article. Nat Protoc 20 (1): 180–219. https://doi.org/10.1038/s41596-024-01045-4.
Luecken, M. D., and F. J. Theis. 2019. “Current Best Practices in Single-Cell RNA-Seq Analysis: A Tutorial.” Journal Article. Mol Syst Biol 15 (6): e8746. https://doi.org/10.15252/msb.20188746.
Saelens, W., R. Cannoodt, H. Todorov, and Y. Saeys. 2019. “A Comparison of Single-Cell Trajectory Inference Methods.” Journal Article. Nat Biotechnol 37 (5): 547–54. https://doi.org/10.1038/s41587-019-0071-9.
Wolf, F. A., P. Angerer, and F. J. Theis. 2018. “SCANPY: Large-Scale Single-Cell Gene Expression Data Analysis.” Journal Article. Genome Biol 19 (1): 15. https://doi.org/10.1186/s13059-017-1382-0.